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Abstract 

Precise centroid estimation plays a critical role in accurate astrometry using telescope images. Con- 
ventional centroid estimation fits a template point spread function (PSF) to the image data. Because 
the PSF is typically not known to high accuracy due to wavefront aberrations and uncertainties in op- 
tical system, a simple Gaussian function is commonly used. PSF knowledge error leads to systematic 
errors in the conventional centroid estimation. In this paper, we present an accurate centroid estimation 
algorithm by reconstructing the PSF from well sampled (above Nyquist frequency) pixelated images. 
In the limit of an ideal focal plane array whose pixels have identical response function (no inter-pixel 
variation), this method can estimate centroid displacement between two 32x32 images to sub-micropixel 
accuracy. Inter-pixel response variations exist in real detectors, e.g. CCDs, which we can calibrate by 
measuring the pixel response of each pixel in Fourier space. The Fourier transforms of the inter-pixel 
variations of pixel response functions can be conveniently expressed in terms of powers of spatial wave 
numbers using their Taylor series expansions. Calibrating up to the third order terms in this expansion, 
we show that our centroid displacement estimation is accurate to a few micro-pixels using simulated data. 
This algorithm is applicable to the new proposed mission concept Nearest Earth Astrometry Telescope 
(NEAT) to achieve mirco-arcsecond accuracy in relative astrometry for detecting terrestrial exoplanets. 
This technology is also applicable to high precision photometry missions. 



1 Introduction 

Consider an Nx N detector array with coordinate (to, n) labeling the pixel in m-th row and n-th column. 
Given the intensities I mn measured by pixels (to, n), to, n = 1, 2, • • • , N, a straight forward estimate of 
centroid is 

x c — ^ J i Vc — ^ J ) (J-J 

where {x mn , y mn ) are the x and y coordinates of the center of pixel (to, n) and the summation is over all 
the pixels relevant to the centroid estimation. Typically only a small array of pixels is used for estimation 
because enlarging the size of array degrades the signal to noise ratio (SNR) significantly. The pixels at 
large distances from the center of the image detect very little photons but are heavily weighted by their 
large coordinates. Estimation ([!]) suffers systematic errors from using a relatively small array of pixels. 
The point spread function (PSF) fitting algorithms pQ supersede this straight forward centroid estimation 
because it avoids amplification of noise from multiplying large coordinates. The main challenge for the 
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PSF fitting approach is the knowledge of the PSF. Computing PSF using a diffraction model requires 
both knowledge of the optical system and wave front aberrations that both are usually hard to obtain. In 
the past, Gaussian functions has been popularly usedpQ However, in order to achieve a highly accurate 
centroid estimation, a more precise PSF is needed. Mighell has done work on PSF fitting using the digital 
images by using 21-point damped sine interpolation [2]. In this paper, we work along the same line to 
reconstruct the PSF from the pixelated images using the fact that the PSF is a bandwidth limited signal. 
Theoretically, bandwidth limited signal can be reconstructed to any accuracy as far as the sampling is 
above the Nyquist frequency and the number of the sample is sufficiently large. Using 32x32 images, 
the truncation error causes less than a micro-pixel error in centroid displacement estimation. 

Due to the complex micro-structure of detectors and charge diffusion effect between pixels, the pixel 
detection response varies over the physical detection area of a pixel[3], which is referred as inira-pixel 
variation. The pixel counts is a convolution of the pixel response function with the photon energy flux 
function. If each pixel has identical response function, pixel counts still represents a bandwidth limited 
function sampled at the pixel grid and thus the sampling theorem is still applicable for reconstructing 
an effective PSF for the detector. Past measurements found dominant intra-pixel detection variation is 
indeed common to pixels [Ij. For micro-pixel centroid estimation, we still need to take into account the 
inter-pixel differences of the pixel response functions. The inter-pixel response variations make the pixel 
counts no longer represent sample values of a strictly bandwidth limited function. To characterize inter- 
pixel variations, we use laser metrology to measure the pixel responses in Fourier space. It is convenient 
to parametrize the Fourier transforms of the pixel response functions in terms of powers of spatial wave 
numbers using their Taylor series expansions. The leading order effect of inter-pixel variations is the 
average pixel response or flat-field response, which can be measured as response to a uniform E-field. 
The first order correction is an effective geometric pixel location shift for each pixel. The second and 
third order corrections are quadratic and cubic polynomials of the spatial wave numbers. As we include 
more terms in the expansion, the model becomes more accurate and the centroid estimation becomes 
more accurate. 

In this paper, we present results based on simulations to demonstrate the capability of calibrating 
the inter-pixel variations for performing micro-pixel level centroid estimations by including up to third 
order terms in the Taylor series expansion of the Fourier transforms of the pixel response functions. This 
algorithm can be utilized by the proposed mission NEAT [7] to perform micro- arcsecond level relative 
astrometry with a one meter telescope and thus detect terrestrial exo-planets in habitable zone. 

In addition to precise astrometry, our pixel calibration technique is also applicable to high precision 
photometry, which requires to characterize the pixel response functions. With the calibration of pixel 
responses and the reconstruction of PSF, the photometry is no longer sensitive to pointing errors and 
detector response variations, which limit the performance. Because we measure the pixel responses in 
Fourier space, it is especially convenient for studying photometry when the images are not too under- 
sampled because we only need to measure the pixel response to a modest spatial frequency range to 
cover the bandwidth of the PSF. 

2 Model and Algorithm description 
2.1 Pixel intensity model description 

To simplify the formulation, we assume that the image is stable over a sampling period so that the 
temporal integration by the detector is simply an overall factor of the duration of the sampling period. 
The image is represented by photo-electron counts recorded by the pixels. A model describing the 
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photo-electron counts recorded by pixel in the m-th row and n-th column is expressed as 



OO 



— OO 



Imn(x c ,yc)= dx dyI(x-x c ,y-y c )Q mn (x,y) (2) 



where Q m n(x,y) is the pixel response function (PRF) of pixel (m,n) to a point illumination at (x,y) in 
the detector planejUE], I(x,y) is an input intensity function, and (x c ,y c ) is the location of the centroid 
of the image. For a point source, which we will consider exclusively, I(x, y) is the point spread function 
(PSF) related to the focal plane E-field via 

I(x,y) = \E(x,y)\ 2 . (3) 

By Fourier optics, E(x, y) is related to the E-field at the input pupil of the telescope Ei(x, y) via a Fourier 
transform 

E(x,y)=Mj dx' J dy'P(x',y')Ei(x',y')ex.p^ — (xx' + yy')j (4) 

where J\f is a normalization factor, A is the wavelength of the light, / is the focal length of the telescope, 
and P(x, y) is the aperture function, whose value is 1 inside the aperture and outside the aperture. 
We will focus on monochromatic light case and briefly discuss about the polychromatic case in section [5] 
Polychromatic case requires one extra integral over the wave number, which makes the formulation 
slightly complicated. The main steps of derivation remain the same. 

Performing a change of variable (k x , k y ) = 2tt/ (A/) (x f , y'), expression Q becomes a two dimensional 
Fourier integral for E(x,y) as 

e( X , y) = m £ * £ ik yP (M t „ M.^) El M.^ ( 5) 

where (k x ,k y ) represents the spatial frequency. Because P(x,y) vanishes outside the aperture of the 
telescope, E(x, y) is a two dimensional bandwidth limited signal. By the theorem of convolution, the 
Fourier transform of I(x, y) is the Fourier transforms of E(x, y) convolved with its complex conjugate. 
Therefore, I(x, y) is also a bandwidth limited function with bandwidth being twice of that of E(x, y) 
from the process of convolution. Let D be the diameter of the aperture of the telescope, the bandwidth 
of I(x,y) is then limited by \k x \ < 27rD/(A/), \k y \ < 2TrD/(Xf)r\ Assuming the pixel size is smaller 
than the Nyquist sampling spacing f\/(2D), by sampling theorem, the intensity function I(x,y) can 
be precisely reconstructed from a pixelated image given an infinitely large detection array. In practice, 
this is complicated by two things. First of all, we can not have infinite samples. Therefore, we have 
to truncate the reconstruction process at some finite size, which introduces truncation error. We found 
that for 32x32 array size (including about the 7th Airy ring) the truncation errors are small enough 
for micro-pixel centroid estimation. Secondly, the pixel intensity model involves pixel response functions 
Qmn(x,y), which depend on pixels. If Q m n(x,y) does not depend on pixels, as we will see later, the 
intensities measured still correspond to sampled values of a new bandwidth limited function at the 
pixel grid. Treating the new bandwidth limited function as an effective PSF, we can reconstruct it by 
sampling theorem. In fact, experiments found that the leading order intra-pixel variation is common for 
all the pixels For micro-pixel level centrioding, however, we can not ignore the inter-pixel variations 
of Q mn(x,y). The pixel counts no longer strictly correspond to sampled values of a bandwidth limited 
function. Fortunately, the inter-pixel variation is only a small fraction of the total response, whose 
Fourier transform can be characterized by laser metrology fringe measurements. 



1 In fact, it can be shown that the non-zero Fourier frequencies satisfies k% + k'i < (2irD /(Xf)) 
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To illustrate this, we rewrite the pixel intensity model as 

coo poo roo 



Imn(x c ,y c )= I dx dy dkj dk y Z(k x ,k y )e i{ - k ^ x - x ^ +k y { y-y^Q mn {x,y) 

D J — OO J — OO J — OO 

/oo 
dk y l(k x , ky)e i ^ kx ^ m+1 / 2 ) a - x o)+k y ((n+i / 2 ) a -yc)) 
-OO 

/oo roo 
^ / dyQ mn ((m+l/2)a+x, (n+l/2)o+y) e^ x+k ^ 
-OO J —OO 
[•OO 

rlk I dk T(k k \C) (k k \p i [ k x{(m+l/2)a-x c )+k y {{n+l/2)a-y c )\ fn\ 



where I(k x , k y ) is the Fourier transform of I(x, y) and Q m n{k x , k y ) is defined by 

/oo roo 
dx / dy Q mn ((m+l/2)o+x, (n+l/2)a+y) e *(fc*+M) . (7) 
-00 J —00 

The pixel response function related effects are all captured by the Fourier transform Q ma {k x , k y ). 
For the case that all pixels have identical response function, 

Q mn ((m+l/2)a+x, (n+l/2)o+y) = Qoo (o/2+ar, a/2+y) ; (8) 

and Qmn(k x , k y ) no longer depends on indices m, n, 

Qmn{k x , ky) = Qoo{k x , ky) . (9) 

Putting this in the expression ([6]), we obtain 

roo roo 

Imn(x c , y c ) = dk x dk y l(k x , k y )Q 00 (k x , ky y{h,{{rn + l/2)a- Xc ) + k y {{n + l/2)a-y c )) (1Q) 



Defining an effective PSF I(x,y) 

roo roo 

I(x,y)= dk x dk y l(k x ,k y )Qoo(k x ,k y )e l ^ x+k yy\ (11) 



—00 j —00 



ImnixcVc) corresponds to values of function I(x,y) at grid points (x,y)= ((m+l/2)a— x c , (n+l/2)a— y c ), 
m,n = 0, ±1,±2, Because I(k x ,k y ) vanishes for \k x \ > 2irD/(Xf) or \k y \ > 2irD/(\f), I(x,y) is 
bandwidth limited and the sampling theorem is applicable for reconstructing the effective PSF I(x, y) 
using I mn (x c ,y c ). 

In general, response function Q mn (x,y) depends on (m,n). It is useful to define an average response 
function Q(x,y) 

Q(x, = Qmn({m + l/2)o + x, (n + l/2)a + y) . (12) 

m,n 

and its Fourier transform 

Q(k x , k y ) ^JdxJ dyQ(x, j,)e*(*-* + M) . (13) 

Because the common intra-pixel response function does not affect the property of being sampled values 
of a bandwidth limited signal, it is convenient to factor out the common pixel response to parametrize 
the pixel dependent terms as 



(1 (h h \ID(h k \ — n p ik x l±x mn +ik y i\ymn i i i ft lc J ' 4- n, U h 

^cmny^xi ">y } / ^yr^xi ">y ) — Hmn^ L \ ' JL mn rx ' X \ Pmn^y \ !mn rh x l ^y 

~\~ Q"mnk x ~\~ b mn k x ky -\- c mn k x ky + d mn ky + ■ • • , (14) 
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where q mn represents the flat field response of pixel (m, n), Ax mn and Ay mn represent effective geometric 
pixel location shifts along x and y directions for pixel (m, n) to deviate from a regular pixel grid location 
((m + l/2)a, (n + l/2)a). a mn , /3 mn , 7 mn are parameters specifying quadratic behavior in the amplitude 
of Qmn(k x , k y ) that are pixel dependent. a mn , b mn , c mn , d mn are third order term coefficients. This 
parametrization is based on a Tayler series expansion of Qmn(k x , k y ) in terms of polynomials of k x and 



2.2 Centroid displacement estimation algorithm 

We now turn to present the algorithm for estimating the centroid displacement between two images. The 
main idea is to resample one of the two images at a grid shifted from the default grid by an offset. When 
the resampled image best matches the second image, the offset at which we resampled the first image 
is the estimated centroid displacement. We call the first image (the image to be resampled) reference 
image. 

Based on pixel grids, we use the following discrete frequencies 

(k 3 , h) = — (j, l),j = -N/2, ■ ■ ■ , 0, 1, N/2 - 1 , I = -N/2, ■ ■ ■ , 0, 1, N/2 - 1 . (15) 

With this, our model is expressed as 

N/2-1 N/2-1 ^ it, u\ 
Imn(Xc,y C )= £ J2 Mj rT { j ' e %l(^+l/2)a- a;e ]+ifc i [(n+l/2)a-y e ] 
j=-N/2l=-N/2 Q(kj,h) 

where 

M jl =X(k 3 ,k l )Q(k J ,k l ). (17) 

Pixel response calibration measures Q m n(k x ,k y )/Q(k x ,k y ) at a set of spatial frequencies covering the 
bandwidth of PSF. It is convenient to parametrize Qmn(k x ,k y )/Q(k x ,k y ) using low order terms in ex- 
pression (14). If keeping the terms up to the third power, we have then parameters q mn , Ax mn , Ay mn , 

dmn- They can be estimated by fitting expansion \\M to calibration mea- 
surements, which are values at a set of spatial frequencies pre-determined by the metrology system for 
pixel response calibration. 

To reconstruct the PSF, we need to estimate Mji- Let us choose the coordinate so that the centroid 
of the reference image is at the origin, i.e. x c = 0, y c = and its pixel intensities are represented by I mn - 
Detector calibration measures flat pixel responses q mn , effective pixel offsets Ax mn , Ay mn , quadratic 
amplitude coefficients a mn , f3 mn and 7 mn , and the third order phase parameters a mn ,b mn ,c mn ,d mn . 
Relation (16) 

N/2-1 N/2-1 ^ (U U\ 

j0 _ Mjl K U c ik i (m+l/2)a+ik l (n+l/2)a 

j=-N/2l=-N/2 Q{kj,ki) 

can be inverted to solve for Mji, which is the frequency representation of the image. We use expres- 
sion (16) to resample the reference image. The centroid displacement of the second image Imn relative 
to the first image can be estimated by solving a weighted nonlinear least-squares fitting as 

(x c ,y c ) = min£ \I mn - I mn (x c , y c )\ 2 ' W mn (19) 

mn 

where W mn is weight factor for pixel (m,n). The choice of W m n affects sensitivity to noise. Here we 
simply use equal weights and save optimizing the weights to achieve best noise sensitivity as a future 
topic. 
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When we use discrete frequencies (15), the signal is parametrized as periodical function in space with 
period Na. For N = 32, the image encircles the 7th Airy ring; the PSF is small enough at the boundary 
so that the truncation error is negligible. 



3 Pixel response calibration using laser metrology 

We measure the pixel response functions in Fourier space by observing the response of the detector to a 
sinusoidal intensity illumination pattern. The sinusoidal pattern is generated by fringes from interfering 
two laser metrology beams far away so that the wavefront is close to a plane wave over the spatial exten- 
sion of the detector. To make the identification of the laser fringes easier, we use Acoustic Modulation 
Oscillator (AMO) to offset one of the laser's frequency from the other by a few Hz so that the fringes 
move across the CCD. Therefore, the illuminating intensity may be modeled as 

I(x, y, t) = i] nct + /£ lct + 2^/ 1 met 7 2 met Ree i(fc " cta ' +fc r t s/+ A ^*) , (20) 

where /f et and If 16 * are the intensities of the two lasers, (&™ et , fc™ et ) is the spatial wave number of the 
laser fringes, and Aui is the angular frequency difference between the two lasers introduced by the AMO. 
The output counts of pixel (m, n) is then 

Cn t (*)= / dx J dyI(x,y,t)Q mn (x,y) 

= (ir et +I 2 mCt )<3mn(0, 0)+2y / 7pi^Re{Q mn (fcr t , ^ et ) e i ( fc S 1Ct ( m + 1 /2)a+fc- t (n+l/2)a+A^)| ^ 

where we have used the definition ([7]) . The temporal variation of the intensity at each pixel is a sinusoidal 
function plus a constant. By estimating the amplitude and phase of the sinusoidal temporal variation, we 
get the complex Fourier transform Q mn (/c™ et , k™ et ) because the overall intensities of the two metrology 
beams If iet and I™ * can be easily measured. By having different separations between the laser beams and 
distances between the laser source and the CCD, we can generate fringes of different spatial frequencies 
and thus measure Qmn{k™ et , k™ et ) at various values of (k™ et , /c™ et ). Note that we only need to measure 
the Fourier transform of the the pixel response functions to the highest frequency in the PSF, which is 
limited by the aperture size of the optics. For Nyquist sampled images, it is only necessary to measure 
the spatial frequency to n/a, where a is the spacing between pixels. 



4 Results using simulated data 



We use simulated data to validate our concept of the algorithm. In subsection 4.1, we first show that 
for an ideal detector whose pixels have the same response function, we can achieve sub-micro-pixel level 
accuracy in centroid estimation. 



4.1 Results for an ideal detector 

In our simulation, the diameter of the telescope D=lm with focal length / = 40m to be consistent with 
NEAT[7]. We focus on monochromatic source with the wavelength A = 600nm. In section [5j we discuss 
the case of polychromatic source. We only consider an array of 32x32 pixels, i.e. the dimension N=32 
because this is sufficient for achieving micro-pixel centroid estimation. The corresponding \/D at the 
focal plane is 24/mi. The pixel size is 10/im., which samples above the Nyquist frequency (l/(12yum)). 
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Pixel response function Pupil phase aberration (rad) used in simulation, RMS=V20 




Figure 1: A common pixel response function for all the pixels used in the simulation (left) and low order 
phase aberration (right) in radian used in simulation with RMS ~ A/20. 



All the pixels have the same pixel response function, which is displayed in the left plot in Fig. [TJ It is 
modeled as a Gaussian function multiplied by low order polynomials (up to 4th order), 



Q(x, y) = exp(-(x 2 + y 2 )/r 2 g ) c + c\x + c 2 y + c 3 x 2 + c^y 2 + c 5 xy + 



(22) 



where x, y are coordinates in the detector plane in unit of pixel, r g = 0.5, and the coefficients Cj are drawn 
from Gaussian random number generators with standard deviation being 0.05. The mean values are all 
except for cq whose mean is 0.8. This pixel response function is roughly similar to the the intra-pixel 
variation for a backside-illuminated CCD [6]. We include A/20 RMS low order wavefront aberrations 
parametrized by the first 15 Zernike polynomials, whose amplitudes are randomly generated. The right 
plot in Fig. [T] displays the phase aberration over the telescope pupil used in our simulation. The two 
plots (left for X centroid and right for Y centroid) in Fig. [2] displays the centroid estimation errors for 
a grid of X and Y offsets between the two images within range [-0.5, 0.5] pixel. These errors are due to 
truncation and are no more than 0.1 micro-pixel. 

4.2 Pixel response calibration and results 

For realistic detectors, pixel response functions vary from pixel to pixel, which we call inter-pixel vari- 
ations. Pixel response calibration measures the pixel response functions of all the pixels. The Fourier 
transforms Q mn (,k x , k y ) of the pixel response functions Q mn (x,y) of all pixels are measured at various 
spatial frequencies. To do this, we illuminate the detector using a cosine intensity pattern. For NEAT, 
this is achieved by interfering two metrology lasers at various spatial frequencies . The fringe pattern 
gives a sinusoidal illumination on the pixel array and the separation of the two lasers determines the 
wavelengths of the fringes. AMO is used to offset the frequency of one laser from the other by a few Hz. 
The temporal variation of the pixel intensities are then used to estimate Fourier transforms Q mn (kx, k y ) 
as discussed in section [3J We use expansion (14) to parametrize Q mn (kx,ky). For micro-arcsecond 



accuracy, we nominally keep the terms up to second order or third order in k x and k 



v 



To simulate the inter-pixel response variations, we make the low order coefficients c, in Eq. (22) 



pixel dependent by replacing a with a + cf 111 , where cf 1 " are random numbers drawn from zero mean 
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X Centroid error, RMS=0.019(fipixel), identical pixel response,W20 




Y Centroid error, RMS=0.0055(^pixel),identical pixel response,^/20 
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Figure 2: Centroid offset estimation errors X (left) and Y (right) in micro-pixel unit for an ideal CCD and 
wave front with A/20 RMS low order aberrations parametrized by the first 15 Zernike polynomials. The pixel 
size is 10/zm. 



Gaussian random number generators with standard deviation being 0.01. We also add 2% white noise 
to the pixel response as a multiplicative factor as well as random geometric pixel location shifts for 
all the pixels along both x and y directions with 0.01 pixel RMS. The left plot in Fig. [3] displays a 
typical intra pixel responses of a 10x10 pixel array at the upper left corner of the 32x32 array. Here 
the interpixel cross talk from diffusion is not displayed to avoid overlap. The zeroth order calibration 
measures the pixel response to a flat field, i.e. Fourier transform of the pixel response function at spatial 
wave number (k x ,k y ) = (0,0). The right plot in Fig. [3] displays a typical flat field calibration results. 
Fig. [4] shows systematic centroid displacement estimation errors with only flat field response calibration 
for different centroid offsets between the two images, whose range is [-0.5, 0.5] pixel along both x and y 
directions. The next level of calibration measures the effective pixel locations deviating from a regular 
grid. The effective location of pixel (m, n) is estimated by fitting Ax mn k™ et + A.y mn k™ et to the phase 
of the estimated Fourier transform Q mn (/c™ et , /c™ et ) at different values of spatial frequency (£;™ et , fc™ et ) 
available from the metrology system. Fig. [5] displays the estimated the effective pixel location deviating 
from a regular grid. Now the errors shown in Fig. [6] are significantly reduced to tens of micro-pixels 
after including the effective pixel locations in estimation. However, the errors are still large for micro- 
pixel level astrometry. We further include the second order amplitude terms in expansion (14) and 
estimate coefficients a 

Q 



i(£:™ et , fc™ et ) at different values of (kf et ,k^ 



Prnn, and 7 mn by fitting expansion (14) to the estimated Fourier transforms 



, ,. . ,„ ,. The second order coefficients are displayed in Fig. |7| 

The two plots in Fig. [8] displays the centroid estimation errors for performing pixel response calibration 
that estimates the flat field response, the effective pixel locations, and the second order amplitude 



corrections. Including the third order terms in expansion (14) in our pixel response calibration enables 



us to achieve centroid estimation accuracy to be a few micro-pixels over [0.5, 0.5] pixel range along both 
x and y directions. The corresponding coefficients are displayed in Fig. [9} Fig. 10 shows systematic 
centroid displacement estimation errors for including through the third order phase terms. The RMS is 
only a few micro-pixels. 
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Figure 3: Simulated intrapixel detection variation for a 10x10 array (left) and flat field response calibration 
result 32x32 (right). 




Figure 4: Centroid offset estimation errors X (left) and Y (right) in micro-pixel unit. The wavefront error is 
A/20 specified by the first 15 Zernikes polynomials. Only a flat field calibration is applied. 
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Figure 5: Effective pixel offsets along x (left) and y directions, the units are pixel 




Figure 6: Centroid offset estimation errors X (left) and Y (right) in micro-pixel unit. The wavefront error 
is A/20 specified by the first 15 Zernikes polynomials. A flat field and pixel effective location calibration is 
applied. 
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Figure 7: Second order amplitudes a mn (left), (3 mn (middle), and j mn (right), coefficients for terms k%,ky, 
and kxky. 




Figure 8: Centroid offset estimation errors X (left) and Y (right) in micro-pixel unit. The wavefront error is 
A/20 specified by the first 15 Zernikes polynomials. Calibration includes second order terms. 
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Figure 9: Third order term coefficients, ak\ (top, left), hk 2 x k y (top, right), ck x k 2 (bottom, left), and dky 
(bottom, right). 




Figure 10: Centroid offset estimation errors X (left) and Y (right) in micro-pixel unit. The wavefront error 
is A/20 specified by the first 15 Zernikes polynomials. 



12 



Allan deviation of the estimated centroid vs number photons 
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Figure 11: Centroid estimation error Allan deviations vs number of photons integrated. 
4.3 Noise sensitivity 

So far, we have not included any noise. We now study the sensitivity to photon shot noise, which is the 
dominant source of random errors. We set the total number of photons for each image to be 10 8 and use 
Poisson random number generator to simulate the photon shot noise. We simulated 1000 images with 
their centroid positions randomly distributed within ± 0.5 pixel relative to a reference image along both 



x and y directions. The reference image have the same level of photon and shot noise. Fig. 11 displays 
the Allan deviation of the estimated centroid displacements relative to a reference image as function of 
the total number of photons integrated. The green dash line shows an empirical sensitivity formula 



a « 1.5/^iV ph (23) 

for uncertainty of centroid estimation using an equal weight in the least-squares fitting, where iV p h is the 
total number of photons. To reach 10 mirco-pixel precision, we shall need about 2.2xl0 10 photons. The 
Allan deviation shows the variation between 1000 images with respect to the same reference image. How 
does the noise in the reference image affect the results? It turns out that the noise in the reference image 
causes mostly an overall offset to all the centroid displacements, or it is an excellent approximation that 

d(A, R{) - d(B, R{) « d(A, R 2 ) - d(B, R 2 ) , (24) 

where d(A,Ri) represents the displacement of centroid of image A relative to reference image R± and 
so on. Fig. [12] shows the centroid displacements relative to two reference images, which are the same 
image with two different realization of photon noise. It is easy to see that the difference of the centroid 
displacements is insensitive to the uncertainties in the reference image ( less than 1 micro-pixel). We are 
performing an equal weight least-squares fitting, it is possible to optimize weights W mn to achieve the 
best sensitivity to noise, which will be a subject for future study. 

4.4 Wave front stability 

In this section, we study the sensitivity of our algorithm to wave front changes, i.e. the two images 



are taken with slightly different wavefront aberrations. The left plot in Fig. 13 shows a low order small 



wavefront change with A/1350 RMS. (High order wave front aberrations generates speckles in the image 



13 



Relative centroid depends weakly on the reference image 
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Figure 12: X centroid displacement error using reference image with different noise realizations. 



plane while the low order wavefront change generates a change in the shape of the PSF. For PSF fitting, 
the centriod estimation is more sensitive to low order wavefront aberrations.) The main effect of a 
wavefront change between the two images is an overall constant, which is not a problem for differential 
astrometry. However, this overall constant depends on the pixel response functions. We shall need to 
show that this pixel dependency in the overall constant is small so that it is common for different portions 
of the detector. To do this, we instantiate a second detector whose pixel response function is slightly 
different from our nominal detector as a result of the random feature of our simulation. 



The right plot in Fig. 13 shows the difference of the pixel responses between the two cameras which 
comes from the randomly instantiated low order polynomials. We now consider the case where the wave- 
front changes by an amount of A/1350 RMS, as shown in the left plot in Fig. 13, between taking the two 
images whose centroid displacement we are interested to estimate. Fig. [14] displays the centroid displace- 
ment errors due to the A/1350 RMS low order wavefront change for a grid of centroid offsets between 
the two images. The main effect of a wavefront change between the two images is an overall constant 
independent of the actual displacement between the two images; the variation is only a few micro-pixels. 
If the overall constant is common for both stars, this effect cancels for differential astrometry. 

Because the images of the two stars are at different locations in the focal plane array, we examine 
whether the different pixel response functions at different portions of the detector coupled with the 
wavefront change leads to significant error. Fig. 15 displays the centroid errors for using the second 
detector. To avoid common error cancellation, we put the reference image at [0.25, 0.25]pixel instead 
of [0, 0] pixel as for the first detector. Again, we can see that the dominant effect of wavefront change 
between the two images is an overall offset in the centroid estimation. Because the overall offset is 
not sensitive to the difference between two detectors, the overall differential displacement caused by the 
wavefront difference cancels leaving a residual of a few micro-pixel RMS. See Fig. [IB) 
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Phase aberration change used in simulation, RMS=V1348 Intra pixel response difference 
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Figure 13: Pixel response difference between two detectors in a 10x10 array (left). Wavefront variations 
between two observations, RMS ~ A/1350, simulated by low order polynomials with randomly generated 
coefficients. 



X Centroid error, mean±std=935.9±1.2(|Apixel) 



Y Centroid error, mean±std=3767±3.2(jipixel) 
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Figure 14: Centroid displacement estimation errors X (left) and Y (right) in micro-pixel unit. Both images are 
taken with a common low order wavefront error of A/20 RMS specified by the first 15 Zernikes polynomials. 
The low order wavefront difference between the two images has RMS ~ A/1350. 
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X Centroid error, mean±std=936.7±0.66(jipixel) Y Centroid error, mean±std=3770±0.72(^pixel) 




X Centroid Offset (pixel) X Centroid Offset (pixel) 

Figure 15: Centroid estimation errors X (left) and Y (right) in micro-pixel unit for the second detector. The 
images are taken with a common low order wavefront error of A/20 RMS specified by the first 15 Zernikes 
polynomials and a differential low order wavefront error of RMS A/1350. 




Figure 16: The difference between the centroid displacement estimation errors X (left) and Y (right) in micro- 
pixel unit for two different simulated detectors. For each detector, the images are taken with a common low 
order wavefront error of A/20 RMS specified by the first 15 Zernikes polynomials and a differential low order 
wavefront error of RMS A/ 1350. 
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X centroid error, RMS=0.92|Apixel, field 1 Pixel size =8.28um Y centroid error, RMS=0.42|Apixel, field 1 Pixel size =8.28um 
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Figure 17: The estimation errors of the centroid offsets for broadband images. 

5 Polychromatic effect 

In this section, we discuss the polychromatic effect. For white light source like stars, image intensity 
model ^ needs an extra integration over the photon wavelengths weighed with the source spectrum, 

/oo poo poo 

dkj dk y / dA5(A)X(^,^,A)Q mn (^,^,A)e t ^« m+1 / 2 ) a -^) +fc ^(" +1 / 2 ) a -^)] 
-oo J — oo JO 

(25) 

where S(X) is the source spectral energy density function and we have included the wavelength depen- 
dencies in Fourier transforms of the monochromatic PSF function I(k x ,k y , A) and the pixel detection 
function Q m n(k x , k y , A). As a leading order approximation, we ignore the spectral dependency in the 
pixel detection function. The broadband image model has the same expression as the monochromatic 
model with the following replacement 

POO 

l(k x ,k y )^ dXS(X)l(k x ,k y ,X). (26) 
Jo 

Because the integral over the photon wavelength is a linear operation, the polychromatic signal is still 
a bandwidth limited signal, whose bandwidth is determined by the shortest wavelength. As far as 
the pixelated images are Nyquist sampled, our algorithms works as for the case of monochromatic 
source. Fig. [17] displays the centroid offset estimation errors for simulated chromatic pixelated images 
pairs displaced by various offsets along x and y directions assuming an ideal tophat pixel response 
detection. We note that the error is sub-micro pixel. A slight complication comes from the dependency 
of pixel response on the wavelength, which requires pixel response calibration using metrology at multiple 
wavelengths. The results in reference [6] for backside illuminated CCD shows weak dependency on the 
wavelengths especially at the longer wavelength. Because we only need to calibrate the pixel to pixel 
variations of the pixel response functions, for differential astrometry, only the spectral dependency in 
the inter-pixel variations of the response functions coupled with the star spectral difference could cause 
systematic errors. We expect this effect to be small in general and can be calibrated using laser metrology 
at a few wavelengths. We defer the study of this to a future work. 
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6 Conclusions 



We have presented a systematic frame work for accurate centroid displacement estimation and detector 
characterization. For an ideal detector with all the pixels having the same pixel response function, the 
PSF can be accurately reconstructed without any detector calibration. For realistic detectors whose pixel 
response functions share a dominant common portion and small pixel to pixel variations, it is possible 
to measure the pixel responses functions in Fourier space using laser metrology fringes. Measuring the 
pixel response in Fourier space is especially convenient for well sampled images. Keeping a few low 
order terms (e.g. 3rd order) in the Taylor series expansion of the Fourier transform in wave numbers, 
we can achieve centroid estimations accurate to a few micro-pixels. This enables micro- arcsecond level 
astrometry using telescope images of stars and thus detect earth-like exo-planets. This method is also 
applicable to precise photometry. 
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